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ABSTRACT 

We present calculations of active galactic nucleus winds at ~parsec scales along with the associated obscuration. 

We take into account the pressure of infrared radiation on dust grains and the interaction of X-rays from a central 
black hole with hot and cold plasma. Infrared radiation (IR) is incorporated in radiation-hydrodynamic simulations 
adopting the flux-limited diffusion approximation. We find that in the range of X-ray luminosities L = 0.05-0.6 
L e dd> the Compton-thick part of the flow (aka torus) has an opening angle of approximately 72°-75° regardless 
of the luminosity. At L > 0.1, the outflowing dusty wind provides the obscuration with IR pressure playing a 
major role. The global flow consists of two phases: the cold flow at inclinations 6 > 70° and a hot, ionized wind 
of lower density at lower inclinations. The dynamical pressure of the hot wind is important in shaping the denser 
IR-supported flow. At luminosities ^0.1 L E dd episodes of outflow are followed by extended periods when the wind 
switches to slow accretion. 
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1. INTRODUCTION 

Active galactic nuclei (AGNs) are among the most fascinating 
objects in the universe: various energetic processes interplay 
in a volume of a few cubic parsecs around a supermassive 
black hole (BH); radiation is likely decisive in connecting the 
accretion and feedback modes of AGNs. Dusty, radiation-driven 
flows are very effective in removing large masses of accreting 
gas, returning it back to the galaxy, and naturally limiting the 
growth rate of a supermassive BH by the feedback with the host 
galaxy. Radiation also provides a nonlinear causal connection 
between distant parts of an AGN itself, coupling together vastly 
different physical scales. Despite the pivotal role of radiation 
in the first-principles modeling of AGNs and the significant 
progress in theoretical modeling over the past 10 years, proper 
incorporation of radiation into numerical models still presents a 
significant computational challenge. 

Well-established phenomenology places AGNs into two dis- 
tinct classes that roughly differ by the presence (type I) or the 
absence (type II) of the broad emission lines in the optical 
and UV spectra. According to the AGN unification paradigm 
(Rowan-Robinson 1977; Antonucci & Miller 1985; Urry & 
Padovani 1995), this dichotomy is an artifact of a pure geo- 
metrical origin. In order to be observed as a type II object, an 
AGN must be viewed at high inclination, in extreme cases close 
to the equatorial plane, while in the case of a type I AGN the 
observer’s line of sight is at lower inclination, i.e., closer to 
the symmetry axis. The presence of geometrically thick toroidal 
obscuration is pivotal to the unification scheme of AGNs. 

The discovery of broad permitted emission lines in the 
polarized spectrum of the nearby Seyfert II galaxy NGC 1068 
(Antonucci 1984; Antonucci & Miller 1985) provides evidence 
that a type I nucleus is harbored inside the obscuring ring of 
matter. The observed polarization is likely due to scattering of 
light emitted by the type I core and then being scattered within 
an outflow emanating from the same nucleus. 


Observations of AGNs and quasars reveal that dust in large 
quantities can be present within a few parsecs from the BH, and a 
significant part of X-rays and UV are reprocessed into infrared 
radiation (IR). For example, in luminous quasars almost half 
of the bolometric luminosity can be reprocessed into IR. This 
results in a 1-100 /rm “hump” in the quasar spectral energy 
distribution (Richards et al. 2006). In nearby AGNs, multi- 
temperature dust distributions are observed directly by the Very 
Large Telescope Interferometer. Such mid-infrared observations 
of the prototypical Seyfert II galaxy NGC 1068 (Jaffe et al. 2004) 
and the Circinus galaxy (Tristram et al. 2007) depict a dusty, 
clumpy plasma exposed to external illumination. 

Various models have been proposed to explain the real or the 
apparent thickness of the torus: 

1. A “Warped disk model”: these are global warps in a 
locally geometrically thin accretion disk (i.e., Phinney 
1989; Sanders et al. 1989). 

2. (i) Models of a quasi- static, rotating torus involving 
“clouds” (i.e., Krolik & Begelman 1988; Beckert & Duschl 
2004; Nenkova et al. 2008). Vertical thickness of the torus 
is supported via a large velocity dispersion of clouds, and a 
small-scale magnetic field is required to maintain elasticity 
of clouds during collisions. 

(ii) Models that rely on IR pressure on dust to keep the 
rotating torus vertically thick (Krolik 2007). Most of the IR 
in this model comes from reprocessing of external X-ray 
and UV radiation. 

3. A “wind torus”: a collar of outflowing matter may provide 
obscuration. 

(i) MHD wind: the possibility is that a global magnetic 
field can be involved in the form of an MHD wind (Konigl 
& Kartje 1994; Elitzur & Shlosman 2006) or directly 
supporting a quasi-static torus (Lovelace et al. 1998), or in a 
combination of radiative and magnetic driving (Emmering 
et al. 1992; Everett 2005; Keating et al. 2012). Note that 
a sufficiently strong (of the order of the equipartition one) 
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poloidal magnetic field can also help to solve the problem 
of how the angular momentum is removed from accreting 
matter; however, the origin of such a highly ordered field 
remains unclear. 

(ii) An outflow (or failed wind) driven by IR pressure on 
dust (Dorodnitsyn et al. 2011, 2012). All these models have 
their benefits and drawbacks. 

Models investigating the properties of a self-gravitating, non- 
radiative, dissipationless torus through V-body simulations have 
been made by Bannikova et al. (2012). 

The accreting dusty gas at ^parsec distances from the BH is 
susceptible to self-gravitating instability; since it is exposed 
to intense radiation, strong gravitothermal instabilities can 
develop. 

In the outflowing wind, the nonlinear phase of thermal 
instability is associated with “clumps” (i.e., Elitzur 2008). 
Presumably such clouds are confined by a small-scale magnetic 
field. Whether they are secondary to the flow and behave more 
like “waves,” i.e., temporary structures evolving on dynamical 
or radiation timescales, or real long-lasting “clumps” remains 
unresolved. All models involving clouds need to explain how to 
avoid too much dissipation through the cloud-cloud collisions, 
which would otherwise raise the temperature to a fraction of 
the virial temperature r v i r g . In the present studies, we do not 
consider the effect of the magnetic field (both ordered and 
chaotic). 

It has long been recognized that the geometrical thickness 
of the obscuration presents a significant problem. In order 
to explain why the torus is geometrically thick, models 
have been constructed that rely on gas pressure, magnetic 
forces, and radiation pressure either alone or in combina- 
tion with other mechanisms. If one relies on gas pressure 
to support vertical thickness of the torus, then virial argu- 
ments suggest that the temperature of such gas should be 
of the order of r v i r g = 2.6 x 10 6 Mj/r vc K, where Mj is 
the BH mass in 10 7 M 0 and r pc is the distance in parsecs. 
Such high temperatures are not compatible with the existence 
of dust. 

As mentioned before, the presence of a small-scale mag- 
netic field can provide necessary elasticity to cloud-cloud col- 
lisions. However, this argument is not without a flaw. It is 
natural to assume that turbulence develops. If so, this leads 
to the magnetic field being of the order of the equiparti- 
tion field. A result is that the dissipation of magnetosonic 
waves will raise the electron temperature to be of the order 
of r vir)g . 

The connection between a toroidal obscuration, UV 
absorption, and warm absorbers can be self-consistently 
assessed by global multi-dimensional simulations. The preq- 
uisite for such an endeavor is to include various physical ef- 
fects and couple them to radiation. Additionally, why such 
calculations are currently beyond reach is that high numer- 
ical resolution is needed to resolve the multiphase behavior 
of the outflowing gas. These are the reasons why different 
works including present studies assess smaller aspects of the 
bigger picture. 

The connection between a broad-line region and dusty out- 
flow was considered by Czerny & Hryniewicz (2011). They sug- 
gested that the gas of the dusty accretion disk is weakly coupled 
to the BH due to radiation pressure on dust, resulting in a dusty 
wind that gives the low-ionization part of the broad-line region. 
Bottorff et al. (2000) considered the relation between the ob- 
served properties of UV and X-ray absorbers and those predicted 


from a magnetohydrodynamical self- similar wind model. They 
conclude that in the case of the well-studied Seyfert 1 galaxy 
NGC 5548 the intrinsically clumped absorber is biased toward 
smaller radii in the case of warm absorbers and is likely located 
at higher distances in the case of UV absorption features. 

Other models are based on the possibility that the pressure of 
the infrared radiation on dust may itself be enough to support 
the vertical thickness of the torus. The idea of a static infrared 
support was explored in Krolik (2007), where an approximate 
analytical model for the torus was constructed showing that 
there is enough radiation force to support a static vertically 
thick torus. However, it was argued in Dorodnitsyn et al. (2011, 
hereafter Paper I) that it is very unlikely that the IR- supported 
torus would be static. It was shown that the equilibrium 
between rotational, gravitational, and radiation forces cannot 
be maintained if the temperature in the torus exceeds that 
of T vnj - 3l2(n/l0 5 M 7 r~ 1 ) 1 ' 4 - 981(n/l0 7 M 7 r- l ) l ' 4 K, 
where n is the number density, necessarily resulting in the 
torus being an outflow. In the following paper (Dorodnitsyn 
et al. 2012, hereafter Paper II), we supported this idea with 
radiation-hydrodynamic simulations of the outflowing accretion 
disk wind, driven by the radiation pressure on dust. The 
distribution of IR was calculated using a flux-limited diffusion 
approximation. 

An important simplification was made in Papers I and II: 
only the infrared-driven part of the wind was considered. The 
argument was that this portion of the wind is far enough from the 
BH so that all the conversion of X-rays to IR already happened 
outside the computational domain. That allowed us to consider 
only one group of radiation (the IR). The radiation temperature 
at the boundary was calculated to match the input energetics 
of soft X-rays. In other words, we did not treat explicitly 
the reprocessing of soft X-rays into IR. Correspondingly, we 
were not able to follow the interaction of the X-ray-heated 
hot component of the flow with the cold, dusty IR- supported 
component. 

In this paper, we calculate the wind structure and explicitly 
treat the interaction between X-rays and gas with subsequent 
conversion of X-rays to IR. This allows us to follow in dy- 
namics the transition from hot, photoionized to cold, dusty and 
IR-supported flow. To make calculations numerically tractable, 
we made several major simplifications in treatment of the trans- 
fer of X-rays and their interaction with the gas. For example, 
transfer of X-rays is reduced to simple attenuation with no mul- 
tiple scatterings taken into account, and we approximately cal- 
culate the heating of dust grains by X-rays. The rest of our 
computational radiation-hydrodynamic framework remains un- 
changed from Paper II: the transfer of the infrared radiation and 
its interaction with matter are treated in a flux-limited diffusion 
approximation. The calculations are time-dependent and 2.5D, 
i.e., three dimensions in cylindrical coordinates and assuming 
axial symmetry. 

The organization of this paper is as follows. After this 
introduction, in Section 2 we review basic assumptions of 
our model, including the description of the radiating fluid of 
the AGN dusty wind by means of radiation-hydrodynamic 
equations. Interaction of the cold and hot components of 
the wind with X-rays is described in Section 3, including 
corresponding approximations. In Section 4 we describe the 
numerical setup and initial and boundary conditions (BCs). The 
results are summarized in Section 5, the obscuring properties of 
the wind and the mass loss rates are discussed in Sections 7 and 
8; the results are summarized in the Discussion; the place of our 
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model to a broader perspective of AGN parsec-scale flows and 
the final comments are given in the Conclusions section. 

2. ASSUMPTIONS AND PHYSICAL MODEL 

The flow in which the continuum radiation plays an important 
role is described by the system of radiation-hydrodynamic 
equations. To first order in v/c, these can be formulated in 
the following form (Mihalas & Mihalas 1984): 


D t p + p V • v = 0, 


( 1 ) 


D t \ = 


1 

P 


Vp + g rad - VO, 


( 2 ) 


pD t 



— pV • v — 4 ttxp# + cxeE + n 2 (T — A), (3) 


pD t 



— V • F — Vv : P + 4tcx?B — cxeE , (4) 


temperature 10 2 -10 3 K, is approximately 10-30 times larger 
than that of the electron Thomson opacities (Semenov et al. 
2003). In this paper, when calculating radiation pressure force 
in Equation (5), we take into account only IR pressure on dust. 
Including UV force (in lines and on dust) is left for future 
studies. We adopt a simple procedure for the calculation of g rad : 
if the temperature of the dust T d > T sub , where T sub is the dust 
sublimation temperature (see further in the text), the force is 
assumed to be zero. 

Our description takes into account two types of ra- 
diation: infrared and X-ray. The transfer of X-rays is 
treated in a single-stream approximation. This implies that 
X-rays are simply traced from the point source (i.e., 
corona) adopting e~ Tx attenuation, where Tx is an opti- 
cal depth with respect to X-ray absorption (Dorodnitsyn 
et al. 2008). The distribution of the IR is calculated in a flux- 
limited diffusion approximation. A detailed discussion about the 
validity of Equations (l)-(4) is given in Paper II. We briefly note 
that the closure relation between F and E is obtained adopting 
the diffusion approximation: 


where quantities related to matter are as follows: p, p, and 
e are the material mass density, gas pressure, and gas energy 
density, respectively, and v is the velocity; quantities related to 
radiation are the frequency-integrated moments (for definitions 
see Appendix A): E is the radiation energy density, F is the 
radiation flux, and P is the radiation pressure tensor; xp> Xe 
are the Planck mean and energy mean absorption opacities (in 
cm -1 ), c is the speed of light, B = oT A /Tt is the Planck function, 
a = a c/4 is the Stefan-Boltzmann constant, and T is the gas 
temperature; other notation includes the convective derivative, 
D t = (d/dt) + v- V, and Vv : P denotes the contraction 
( djVi)P l E Note that all the dependent variables in (l)-(4) 
are evaluated in the comoving frame. Note that E = aT 4 , 
where T r is the radiation temperature. Interaction between 
X-rays and the gas is incorporated in the last term of the gas 
energy equation (3), where T and A are, respectively, heating 
and cooling rates. The equation of state for the gas is assumed 
to be polytropic: p = Kpy , where y = 1 + l/n, n is the 
poly trope index, and p = (y — l)e. The gas is one-component, 
one-temperature T = p/x/ plZ, where /x is the mean molecular 
weight per particle, 1Z = 8.31 x 10 7 erg K -1 g -1 is the universal 
gas constant, and plasma with y = 5/3 is assumed to constitute 
the flow. Three components of the flow velocity v R ,v z , and v $ are 
calculated in cylindrical coordinates z, R , assuming azimuthal 
symmetry, 3^ = 0. 

Frequency-independent moments E , F, which appear in the 
above set of radiation-hydrodynamic equations, are obtained 
by calculating angular moments from the frequency-integrated 
specific intensity, /( r, £1, v, t), and are given in Appendix A. 
Forces that are explicitly taken into account include radiation, 
grad pressure, and gas pressure, p -1 Vp. 

These equations should be supplemented by the radiation 
force, g ra d, which is calculated from the following relation: 


F = -DV£, (6) 

where if optical depth is large, r 1, the diffusion coefficient 
reads 

D = cX, (7) 

where X = 1 /(jcp) is the photon mean free path and k = x! P- 
From Equations (5)-(7), one can see that in the diffusion regime 
the radiation force does not explicitly depend on the opacity. 

The diffusion approximation tacitly assumes that optical 
depth x > 1 . Most of the torus where IR pressure is important 
indeed has r d > 1 , where r d is the optical depth of the gas-dust 
mixture in the infrared. 

If r < 1 , the diffusion approximation should be modified so 
a correct limiting behavior at r < 1 is regained. One can see 
that without such modification, when r 1 , the mean free path 
X — > oo, D — >► oo, and |F| — > oo, which is in contradiction 
with a free- streaming limit, where it should be |F| -> cE (7). 
That is, when optical depth becomes small, or when p -> 0, the 
standard diffusion approximation is no longer applicable. 

In order to correctly describe regions of small r, the standard 
approach is to adopt the flux-limited diffusion approximation 
(Alme & Wilson 1974; Minerbo 1978; Levermore & Pomraning 
1981). In this approximation, X is replaced by X* = X A, where 
A is the flux limiter. The flux limiter that we adopted in Papers I 
and II and in the current work is that of Levermore & Pomraning 
(1981): 


6 + 37?lp + R PP 

where Rlp = X|V£|/£. One can see that if r — > 0, then 
R lp -> oo and |F| ~ c E; if r 1, R LP -> 0 and A -> 1/3. 

3. INTERACTION OF RADIATION AND THE WIND 


where xf = Xa + Xt is the total flux mean opacity consisting of 
absorption opacity, Xa> and the Thomson scattering opacity, xt- 
As in Papers I and II, we did not differentiate between xf, Xp> 
and xe- Radiation pressure on dust is most important from the 
dynamical point of view (Paper I). Note that in the infrared, 
the Rosseland mean opacity, /c d = Xd/P of dust with the 


3.1. Photoionization Equilibrium 

The premise of the current work is that the deposition of 
energy at parsec scale is dominated by external sources, i.e., by 
UV and X-rays from the central BH and inner accretion flow. 
Due to the high opacity of the dust-gas mixture, a significant 
part of the UV and soft X-rays is absorbed and reprocessed into 
infrared in a thin layer of thickness, 8l/R\ pc ~ 1.3 x 10 ~ 3 n^ 1 
(Paper I). 
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The dynamical time within the flow is usually much larger 
than the characteristic time of the photoionization and recom- 
bination. Thus, the ionization balance is determined by the 
condition of photoionization equilibrium. The condition of a 
photoionization equilibrium in the wind implies that the state of 
the gas can be reasonably parameterized in terms of one “ion- 
ization parameter,” § — the ratio of radiation energy density to 
baryon density (Tarter et al. 1969). In a popular form, it reads 

% = 4 it FJn ~ 4 x 10 2 • / x T M 6 /(N 23 r pc ), (9) 

where F x is the local X-ray flux, L x is the X-ray luminosity of 
the nucleus, n is the gas number density, N 23 is the column 
density in 10 23 cm -2 , f x is a fraction of the total accretion 
luminosity L B h available in X-rays, and T is a fraction of the 
total Eddington luminosity L E dd = 1.25 x 10 44 M b , where M b 
is the mass of a BH in 10 6 M©. We adopt a simple attenuation 
model in which the local deposition of energy is proportional 
to F x = L x e~ Tx /(4jtr 2 ), where the optical depth at soft X-rays 
is integrated over a straight line connecting the point-source 
corona with a given parcel of gas z x = J K x p dl , and k x is the 
X-ray opacity. 

3.2. Hot and Cold Gas 

To calculate photoionization balance, we take into account 
Compton and photoionization heating and Compton, radiative 
recombination, bremsstrahlung, and line cooling. The rates of 
these processes are calculated making use of the XSTAR pho- 
toionization code (Kallman & Bautista 2001) for the incident 
spectrum, which is a power law with energy index a. As the 
direct incorporation of the XSTAR subroutines into the hy- 
drodynamics code is not feasible, the key ingredient of our 
method is to make use of approximate analytic formulae for 
the heating and cooling rates. These formulae approximate 
extensive tables of heating-cooling rates obtained from XS- 
TAR; these approximations are similar to those of Blondin 
(1994), although modified by Dorodnitsyn et al. (2008) to in- 
corporate better atomic data and power-law ionizing continuum 
with a = 1 : 

r IC (ergcm 3 s" 1 ) = 8.9 x 1(T 36 £ (r x - AT) (10) 
for the Compton heating-cooling, 
r x (erg cm 3 s -1 ) = 1.5 x 1(T 21 f 1/4 T~ t/2 (T X - T)T~ l (11) 
for the photoionization heating-recombination cooling, and 

A(ergcm 3 s -1 ) = 3.3 x 10 _27 T 1//2 

+ (4.6 x 10“ 17 exp(— 1.3 x 10 5 /T) 
x § ( “°- 8 “ a98a) r-l/2+ 10 _24 )5 

( 12 ) 

for the bremsstrahlung and line cooling. Formulae (10)-(12) 
were found to be in reasonable (25%) agreement with numerical 
results. 

A variety of physical and chemical processes influence the 
transformation of X-rays into IR if dense dusty, molecular 
gas is illuminated by external X-rays (Maloney et al. 1996; 
Krolik & Lepp 1989). Rather than attempt to incorporate all 
these processes in our radiation-hydrodynamic framework, in 
the current paper we adopt several approximations that repro- 
duce the qualitative behavior of low-temperature gas exposed 


to external X-ray radiation. At column densities N smaller than 
10 24 cm -2 , photons with energies below 10 keV interact with 
the gas through photoelectric absorption, while Compton scat- 
tering dominates at higher energies. The total ionization per 
particle at a given point in the wind is proportional to § e ff> 
the ionization parameter that takes into account the frequency- 
dependent attenuation of X-ray flux. Maloney et al. (1996) 
showed that heating and cooling rates of such X-ray-illuminated 
molecular gas depend on the modified ionization parameter, 
£ eff = %/N ®2 = 1.26 x 10 ~ l F x /(n^N 22 ), and several regimes 
are approximately distinguished: (1) the regime of highly ion- 
ized gas § e ff » 1. In this regime we make use of approxima- 
tions for the heating and cooling rates (Equations (10)— (12)); 
(2) if 10 -3 § e ff < 1, the gas is primarily atomic, and we 
also adopt heating-cooling rates (Equations (10)— (12)); (3) if 
the ionization parameter is small (§ e ff § m = 10 -3 ), the gas 
is largely molecular. This is where the most complexities due 
to chemistry and radiative interaction with the dust occur. The 
proper incorporation of this regime calls for consideration of all 
the radiative, chemical, and dust network of processes. This is 
challenging even without coupling to hydrodynamics. 

In the regime when § e ff < §m> we assume that the molecular 
gas is in radiative equilibrium with dust, i.e., T g = T d , where 
equilibrium dust temperature is found from the attenuated 
X-ray flux, F x . Thus, we approximately take into account that 
dust can directly reprocess X-rays to IR contributing to the IR 
energy density, E. We approximately take this contribution into 
account, using equation (see Appendix B) 

dT v 

— =n d ccr d (T d - T r ) . (13) 

dt 

Dust temperature is found from the approximate relation T d = 
(4 Fx/(ac )) 1 / 4 (note that E = aT r 4 ), where Fx is the local 
attenuated X-ray flux. If the dust temperature T d > T sub , where 
T sub is the dust sublimation temperature, then it is assumed that 
dust is destroyed and no conversion to IR occurs. 

In our approach, the main difficulty in calculating the energy 
density of IR from X-rays results from a single fluid approx- 
imation. By having only one species (gas) and not evolving a 
separate dust component, we can mimic the interaction between 
dust and X-rays and calculate T d by means of Equation (13), 
while additionally assuming the relation between T d and T in 
a dusty-molecular phase, in which, for simplicity, we assume 
T = T r , with the latter calculated from Equation (13). Other dust 
characteristics in Equation (13) are as follows: n d = pf d /m d is 
the dust grain number density, where /d is the dust-to-gas mass 
ratio and m d is the dust grain mass, and a d = nrj is the dust 
cross section, where rd is the dust grain radius. Throughout our 
calculations, we adopt the following fiducial values for the pa- 
rameters of dust: typical dust grain radius rd = 1 x 10 -4 cm, 
the density of the dust grain p d = lgcm -3 , f d = lx 10 -2 , 
and above the sublimation temperature T sub = 1500 K the dust 
is assumed to be destroyed. Equation (13) for T d (F x ) approx- 
imately describes how X-rays are converted into IR through 
reprocessing by dust. 

If feff ^ and T < c T§^ b , we assume that the exchange 
between T V (E) and T(e) is described by the balance between 
D t (E/p) and the last two terms in Equation (4). The gas 
temperature T can contribute to IR and vice versa. For example, 
gas can contribute to the X-ray-IR conversion through cooling. 
Cooling can be due to radiative processes, i.e., the low energy 
limit of Equations (11) and (12), and due to “adiabatic” cooling 
due to fluid motion. 
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The update of E from Equation (13) is performed adopting 
operator splitting in the source step, i.e., before the update 
of Equation (4) is done. During the update over dt , the new 
E v (t + dt) is calculated from Equation (13). 

4. METHODS 

To solve the system of Equations (l)-(4), we expanded the 
radiation-hydrodynamic framework described in our previous 
papers. The additions include the interaction of X-rays with dust 
grains and with gas. The coupling between the X-ray heating 
and cooling with hydrodynamics has been implemented and 
tested in our earlier works (Dorodnitsyn et al. 2008). Here, we 
have added physics that includes the direct interaction of X-rays 
with dust grains. 

The radiation-hydrodynamic equations are solved using the 
ZEUS code (Paper II) with extensions taken from Dorodnitsyn 
et al. (2008). The radiation-hydrodynamic part is designed and 
tested in Papers I and Paper II and conforms with methods 
and structure of the original ZEUS code (Stone & Norman 
1992). 

In the following, we briefly review the organizational struc- 
ture of the code. Similar to the way pure hydrodynamics is 
treated in ZEUS, the radiation part is also split into the source 
and transport steps. Artificial viscosity is adopted to resolve 
discontinuities in the flow. In the source step the difference 
equations for 3E/3t and de/dt are solved, and in the transport 
step the advection of E and e is calculated. The update of E and e 
is calculated adopting a fully implicit scheme (see Paper II); this 
includes the solution of the diffusion problem and the e E 
update, i.e., the update of gas and radiation energy densities due 
to radiation-gas interaction. 

Similarly, in the X-ray-matter interaction step, heating 
y x -> e and cooling e -> y x due to Compton, radiative re- 
combination, bremsstrahlung, and line cooling are also done 
fully implicitly in a separate update. 

We make use of the R, z cylindrical grid, which extends 
from R m to R om in the radial and from z\ n = 0 to z out = Rout 
in the vertical direction. The lower boundary at Zm = 0 is 
intended to represent a thin accretion disk, which serves as 
an infinite reservoir of gas and dust. Wind can escape from 
all other boundaries, in principle, though the centrifugal barrier 
effectively prevents escape from the inner boundary. A fixed grid 
of size 200 x 200 is adopted for all calculations. Hydrodynamic 
BCs are as follows: at the left (innermost), right, and upper 
boundaries we adopt outflowing BCs. 

At the equator, we fix the inflow speed vo j at a small fraction 
of the speed of sound, leaving the density to adjust according 
to the continuity equation. Hereafter, index j varies over R 
and index i over z. We allow pi s j to vary and fix v is j, that 
is, p isJ = Vi S+ ijPi S+ ij/Vi S and v is = const, provided that 
Vis ^ v is+ i ;S , where v is+ i^ is the sound speed; everywhere in 
the simulations we assumed iWUs+i,s = 0.01. Note that this 
is different from the implementation of the equatorial BC in 
Paper II, where it was assumed that p is j ~ RJ P and velocity 
was found from the continuity equation and only checked after 
the fact that i Vi S ,s- The present implementation of BCs 
allows the disk to be flexible choosing mass loading of the 
wind. The initial distribution of density at the equator was taken 
to be porr 1 , where po = m p n o is the characteristic density 
and nip is proton mass. We find that numerical solution no 
longer resembles the initial equatorial distribution of p, and 
that solutions are not sensitive to the choice of p. Theoretically, 


since p is j is allowed to change, the final solution should only 
weakly depend on initial equatorial distribution of p. Practically, 
however, this is true only for those parts of the disk where the 
wind is developed. We also find that the present implementation 
of BCs produces a more stable solution near the equator, thus 
somewhat relaxing the restriction on the hydrodynamical time 
step. 

The radiative BCs are free- streaming F ~ cE at all bound- 
aries except at the equator, where zero flux BCs are imple- 
mented. 

5. RESULTS 

5.7. Dynamics of the Wind 

Deposition of radiation energy into the flow is controlled 
by the parameter r\ = L bo i/LEdd- We calculate models for 
t] = {0.05, 0.1, 0.3, 0.4, 0.5, 0.6}. The innermost boundary of 
the computational domain is located at radius Ro = 0.5 pc. 
Other parameters are Mbh = 1 x 10 7 Af 0 , fx = 0.5, and 
the initial scale density no = 1 x 10 8 cm -3 . The evolution 
of the radiative flow was followed for 10-30 dynamical times, 
to ~ 1.5 x 10 11 rpl 2 M~ l/2 s. 

In Papers I and II the temperature of the IR was prescribed 
at the boundary in a physically motivated way (to match the 
energetics of X-rays, / x L bo i)- The results obtained in Paper II 
are generally in good agreement with the present modeling. 
However, explicitly including X-rays in calculations reveals 
several important features of the flow that were absent in 
previous studies. 

One of the results of the present study is that the pressure of the 
hot photoionized component is very important in shaping the IR 
flow. Excitation of such a hot, photoionized wind automatically 
mirrors the creation of the IR flow. The interface between the two 
components of the flow is seen in all simulations presented in this 
paper. In all cases, it takes the form of oblique shocks or contact 
discontinuities. A supersonically moving hot gas obliquely hits 
the slowly moving IR flow, creating a discontinuity (in most 
cases an oblique shock wave). 

5.7.7. Model with L = 0.6 L Edd 

Figure 1 shows the development of the radiation-driven disk 
wind for r = 0.6. The outflow develops quickly after about one 
dynamical time, to — 4 x 10 4 yr. The mass-loss rate is (M) — 
0. 1-0.2 M q yr -1 on average and peaks at (M) — 1 M© yr -1 in 
the period of time between 4 x 10 3 yr and 5 x 10 3 yr. The total 
computation spans from 0 to 8 x 10 3 yr. These averaged peak 
numbers include both failed wind and wind that has enough 
energy to escape to infinity. 

The energetics of the wind is measured by computing the 
kinetic output of the wind at the outer boundary of the computa- 
tional domain X: i^n — 2L\^ n /M, where pv 3 / 2 dlL is 

the kinetic luminosity of the wind. For a radiation-driven wind 
we expect low values of L k i n /L bo i; here we have L k i n /L bo i ~ 
5 x 1(T 5 to 2 x 1(T 4 . 

The averaged dynamics of the wind is described by the 
average bulk velocity of the flow (v) = f v pvdV / f y pdV , 
where V is the total or partial volume occupied by the flow. 
Here, the dense part of the flow has (v z ) — a few km s -1 and 
(vr) — a few xlOkms -1 . The maximum velocity reaches 
1000 km s -1 in R and over 600 km s _1 in the z -direction. 

Figure 2 shows a color plot of the distribution of the gas tem- 
perature. The gas spans a vast range of temperatures from being 
in a cold molecular state to a photoionized, high-temperature 
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Figure 1. Color plot of the density, log p, in g cm 3 for the model with L = 0.6 LecM shown at different times given in years. Axes: z: distance from equatorial plane 
in parsecs; R: distance from the BH in parsecs. 

(A color version of this figure is available in the online journal.) 
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Figure 2. Color plot of the gas temperature, log T, in K for the model with L = 0.6 Lsdd* after 4.2 x 10 5 yr. Axes: z: distance from equatorial plane in parsecs; R\ 
distance from the BH in parsecs. 

(A color version of this figure is available in the online journal.) 


warm absorber flow. An apparent feature of our models is that 
the hot wind consists of large-scale inhomogeneities. The cold 
flow does not show such large-scale structure. Note that we 
have a very simple test if the gas is in the cold, molecular-dusty 
phase, § e ff < §m> and that if the hot component is not shown, the 
density structure is much more pronounced in the temperature 
plot. 


It is of interest to address the question of whether there can be 
cold gas in a hot wind as well. However, the limited resolution of 
our studies does not allow us to provide a reliable answer to this 
question. We do not find the coexistence of hot and cold phases 
of gas in any appreciable quantities, but this does not exclude 
the possibility of such coexistence on length scales smaller than 
we can resolve. 
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Figure 3. Color plot of the density, log p, in g cm -3 for the model with L = 0.5 L Edd , shown at different times given in years. Axes: z: distance from equatorial plane 
in parsecs; R: distance from the BH in parsecs. 

(A color version of this figure is available in the online journal.) 


5.7.2. Model with L = 0.5 L Edd : The Importance of 
the Hot Flow Pressure 

Figure 3, which shows a color plot of the density at dif- 
ferent times, is very similar to Figure 1. One can distinguish 
approximately three regions in the flow: (1) a disk-like, geo- 
metrically thin and dense region extending to R ~ 1.5 pc; (2) 
the IR-driven flow with lower density, separated by contact dis- 
continuity from (3) a hot, photoionized flow. As L decreases, so 
does the energetics of the wind. However, the vertical thickness 
of the dense IR-supported flow is determined by the delicate 
balance between (a) the pv 2 pressure of the hot component 
from above and (b) the pressure of IR-supported wind from 
below. 

Decreasing T from 0.6 to 0.5 causes the pressure of the hot 
flow to be less than in the solution shown in Figure 3. The 
balance between the pressure of the photoevaporated wind and 
IR-driven dusty flow is slightly shifted toward IR-dominated 
flow. 

Figure 4 shows the two-dimensional distributions of the 
velocity components v z and vr. The slowly outflowing wind 
is clearly seen; its well-defined boundary extends to 0.5 pc. 
The flow consists of a “core,” i.e., slow-moving wind with 
radial velocities of a few xlO < 100 kms -1 , and v z much 
smaller, of the order of a few 1 kms -1 . The fast lower-density 
component of the flow occupies most of the domain having 
velocities of the order of a few x 100 < 700 kms” 1 . The 
overall characteristics of the velocity field are similar to those 
of the model with L = 0.6 LEdd- Not surprisingly, pumping 
less energy initially in X-rays generates less powerful wind: 
^kin/^bol —lx 10 -5 to7 X 10 -5 . 

The mass-loss rate (M) — 0.1M©yr -1 with the peak 
(M) — 1.1 Mq yr -1 in the period of time between 8 x 10 3 yr 
and 10 4 yr, and the computation spans from 0 to 8 x 10 3 yr. 


5.1.3. Model with L = 0.3 L Edd 

Figure 5 demonstrates the evolution of the density for this 
model at different times. The disk- wind system goes through 
several stages: the domination of the hot wind, puffing up the 
IR-supported disk, squeezing of the latter toward the equator 
by the hot wind, building a high-density disk-like outflow, etc. 
As in previous examples, the system does not come to a quasi- 
stationary state, instead going through such episodes all over 
again. 

An interesting feature that is seen at t = 2.5 x 10 5 yr is in 
fact the residue from the contact discontinuity that existed at 
earlier times. The average velocity (v) — a few x 172 kms -1 . 
This is about 60% of the escape velocity, C/ esc ; the numbers are 
given for t = 20^d yn - The maximum velocity of the dusty flow 
is 242 kms -1 . 

The value of the mass-loss rate is remarkably similar to 
models with higher T: (M) — 0. 1-0.2 M© yr -1 reaching 
1 M q yr -1 . 

5.1.4. Models with L < 0.1 L Edd : Slow 
Accretion with Episodic Outflow 

When the radiation input falls below an approximate thresh- 
old value, no rigorous outflow is observed in our simulations. 
Instead, we see that episodic outbursts of the hot evaporative 
flow are excited from the inner parts of the disk. Here, we found 
this threshold luminosity to be L ~ 0. 1 LEdd - The results for the 
model for L = 0.1 LEdd are shown in Figure 6. 

The density of this wind is too low to provide any considerable 
shielding from X-rays. Without much shielding, it is difficult to 
form a cold and dusty phase that can be accelerated by IR. Such 
an episodic outburst of the dense wind is followed by a gradual 
fallback of the dense component mostly in the z-direction toward 
the equatorial disk. 
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Figure 4. Surface plot of the z- and /^-velocity components. The model is shown with L = 0.5 L^dd and, t = 6.5 x 10 5 yr; R : distance from the BH in parsecs; z: 
distance from the equatorial plane in parsecs. 

(A color version of this figure is available in the online journal.) 
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Figure 5. Color plot of the density, log p, in g cm -3 for the model with L = 0.3 L^dd shown at different times given in years. Axes: z: distance from equatorial plane 
in parsecs; R: distance from the BH in parsecs. 

(A color version of this figure is available in the online journal.) 


The velocity of the hot wind, v < t/ esc , skirts the denser flow 
and falls back toward the equator at larger radii, R ~ 2.5 pc. 
The density piles up there, and the shielding from X-rays 
rapidly increases. Correspondingly, the ionization parameter 
drops below § m , and the cold IR-driven wind again has a chance 


to develop. Note that the vertical component of the radiation 
force is roughly g e ff, z ~ (dE/dr) cos 0, where r is the spherical 
radius and 6 is the inclination angle, measured from the z-axis. 
Since the vertical extent of the IR wind is small, g e ff, z is small, 
and the dynamic pressure of the IR wind cannot overwhelm the 
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Figure 6. Color plot of the density, log p, in g cm -3 for the model with L = 0.1 LecM shown at different times given in years. Axes: z: distance from equatorial plane 
in parsecs; R: distance from the BH in parsecs. 

(A color version of this figure is available in the online journal.) 



pressure of the hot component at higher z- The IR wind stalls 
and starts to slowly fall back (accrete) toward the equator. The 
density outside this dense IR-supported “pancake” drops, and 
the gas gets overionized again. The cold disk itself shrinks from 
both inside and outside. 

We also calculated a model with L = 0.05 L E dd and found 
that the general behavior of this model is similar to the previous 
one with L — 0.1 L E dd- The wind is episodic and intermittent 
with the “IR-supported pancake” phase. 

It is difficult to access the dynamical properties of models 
with L < 0.1L E dd because our equatorial BC excludes the 
influx of matter into the accretion disk. Strictly speaking, by 
assuming the thin accretion disk as a source of matter for the 
wind together with outflowing BC conditions, we assumed that 
the wind, even if very weak, is always present. We will further 
elaborate on this limitation in the following section. 


6. OBSCURING PROPERTIES 

We calculate Thomson optical depth r q = K e pdl, where 
/ is measured along the line of sight at the angle 0 from the 
vertical axis, and averaging is made over all calculated models 
for particular F. Processing different models, we are interested 
in an angle 0 p h where r p h = r(0 p h) = 1, i.e., when the wind 
becomes opaque. 

In the following we summarize results, listing models in the 
order from low to high luminosities. When T = 0.05, the optical 
depth rises sharply from 0.1 to 0.5 at (6) ~ 75° (hereafter 0 e dge)> 
and the photosphere is at (0 p h) — 78°. 

Increasing the luminosity to T = 0.1 causes the dense part of 
the wind to become thicker in vertical extent, optical depth rises 
at the angle 0 e dge — 73°, and the photosphere is at (0 p h) — 75°. 


At T = 0.3 optical depth rises sharply from 0.1 at (0) ~ 73°, 
(0 ph ) — 75°, and 6> e dge — 73°. 

For the model with T = 0.5 the photosphere is at (0 p h) — 75° 
and the torus edge is at 0 e dge — 73°. 

Figure 7 shows a color plot of a radial Thomson optical depth 
measured from a BH to a given point of the computational 
domain. The model shown has L = 0.6L E dd- This is the most 
luminous model we have calculated. It has the photosphere 
located at (0 p h) — 72° and the torus edge at 0 e dge — 71°. 

7. MASS-LOSS RATE 

The opacity of dusty plasma in the infrared domain is 
10-20 times that of electron scattering (Semenov et al. 2003). 
Such huge opacity makes it possible for the AGN to have a 
strong wind, radiating at a fraction of Eddington luminosity 
(defined with respect to Thomson scattering; see Paper I for an 
estimate of an outflow onset). 

Mass-loss rate from the IR-driven portion of the wind can 
be estimated by combining the results from a theory of stellar 
winds and some of the approximations adopted in Paper I. From 
the stellar wind theory, the mass-loss rate from a spherically 
symmetric wind can be found from the following approximate 
relation (Lamers & Cassinelli 1999): 

M/(4tt) ~ L m /c (r IR - l)/(4;r i;°°r IR )r w , (14) 

where T IR = L x /^Edd,iR» v °° is the wind terminal velocity, and 
r w is the wind optical depth in IR. 

For simplicity, as in Paper I, in this section we assume that all 
incident UV and X-ray radiation is reprocessed within a narrow 
conversion layer in which effective temperature, r e ff, can he 
found from the following approximate relation: a Tijx /Add = 
a T e f t , where T = L/L E dd> = 0-5 is the fraction of X-ray 
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Figure 7. Color plot of the radial Thomson optical depth, log r, from a BH toward a given point. Axes: z: distance from equatorial plane in parsecs; R: distance from 
the BH in parsecs. 

(A color version of this figure is available in the online journal.) 


radiation from the total radiation, and the fraction a — 0.5 of 
the incident flux is re-emitted outwardly. Adopting the above, 
one obtains Tir = d //c e , which is approximately 1.25 

for M B h = 10 1 M q , L = 0.5L ddd , and /c d = 10/c e . The 
value of v°° is estimated while neglecting P g in favor of the 
radiation pressure: v°° = y^2 GM/r 0 (r IR — 1) ~ 207kms _1 
adopting the above set of parameters plus the wind launching 
radius Ro = 0.5 pc, and assuming that the wind occupies a 
wedge-like (in a quarter of the domain) space of an opening 
angle of 75°, we arrive at the estimate M ~ 3r w M 0 yr _1 . 
Monte Carlo simulations show that the factor t w can reach 1-5 
for an AGN radiating at a fraction of L Edd (Roth et al. 2012), 
meaning the potential for the wind tori to reach a mass-loss 
rate of 10 M 0 yr _1 . We do not see such mass-loss rates in our 
simulations, which may indicate that instead of being blown 
away, as it certainly would in a spherically symmetric case, 
the gas is effectively ablated toward the equatorial plane by the 
pressure of the hot gas component. As the latter is itself the 
result of interaction with radiation, one can say that radiation 
effectively clears its way through the opening of the torus. 

8. DISCUSSION 

The following average properties of our wind solutions can 
be summarized: 

1. The total mass-loss rate in the IR-driven and warm ab- 
sorber flows is between (M) — 0.1M o yr -1 and (M) — 
1.5 Mq yr -1 depending on the time of the evolution of the 
wind. This also includes the failed wind. The amount of 
gas that has enough energy to escape to infinity is in good 
accord with the energetics derived from accretion. Note 
that approximately 1.6 x 10 _1 tol.3 x 10 -2 M© yr -1 are 
required to be accreted through a geometrically thin accre- 
tion disk to support our models. There is also no problem 
of depletion of the torus since we have an infinite reservoir 
of matter in a razor-thin equatorial disk. 


2. The IR-driven wind has a velocity in the range of a few x 
10-200 km s -1 , and the fast hot wind has a velocity in the 
range of 400-800 km s -1 . 

3. The transmission properties are similar between different 
models. The characteristic angle at which the wind becomes 
opaque to Thomson scattering is about 72°-75° and is 
determined by the balance between the hot, photoionized 
wind and the cold wind supported by the IR pressure on 
dust. It is interesting that the geometry of the torus is close to 
one obtained in Dorodnitsyn et al. (2008) at late times when 
the pressure of the hot, evaporative component becomes 
important. 

4. Both hot and cold components of the wind are 
non- stationary. 

Models such as Konigl & Kartje (1994), Elitzur & Shlosman 
(2006), Emmering et al. (1992), Everett (2005), Keating et al. 
(2012), and Bottorff et al. (2000) are developed on the premise 
of (1) the existence of a strong global magnetic field and (2) 
using the prescription of self- similarity. To a certain extent 
this limits the value of the direct comparison of our results 
with theirs. Indirect comparison would include performing a 
detailed analysis of the spectral properties of the IR-driven 
“wind torus” and performing our simulations for the parameter 
range attributed to a well-studied source (such as NGC 5548, in 
the case of Bottorff et al. 2000). Those are the next things to do 
in order to assess the validity of the current model, and both of 
these goals are a natural continuation of the present studies. We 
will elaborate on observational properties and the implications 
of the time-dependent behavior of our solutions in a following 
paper of this series. 

9. CONCLUSIONS 

Our studies of AGN infrared-driven winds have proceeded in- 
crementally. In Paper I, we considered a toy model that roughly 
approximated an infrared-driven obscuring wind in AGNs. An 
important prediction that resulted was that if the temperature of 
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the dusty plasma in the rotating toms exceeds some character- 
istic value, the IR pressure on dust will inevitably overwhelm 
gravity, creating an outflow. However, this effective tempera- 
ture, which approximately reads as T e ff = ( GMp/ar ) 1 / 4 ~ 
312(ft 5 M 7 /r pc ) 1/4 — 987 (ft 7 M 7 /r pc ) 1/4 K, was derived analyti- 
cally from a simplified model. 

In Paper I, we combined calculations of stationary one- 
dimensional motion in the z -direction with two-dimensional 
calculations of the radiation field. The latter was calculated 
adopting a two-dimensional flux-limited diffusion approxima- 
tion. In Paper II, we relaxed the assumption of one-dimensional 
stationary motion and adopted a 2.5D time-dependent picture. 
One important “toy-model” simplification still remained, how- 
ever; there was no X-ray radiation explicitly taken into account. 
Instead, we calculated the effective temperature of the gas at 
the boundary of the computational domain. In the present work, 
we relaxed this assumption and included X-rays in the global 
picture explicitly. 

Despite these advances, our work still does not consider some 
potentially important physical processes. 

The assumption of a razor-thin disk is idealized. An accretion 
disk at parsec scales is likely self-gravitating, possibly affected 
by gravitothermal instability (i.e., Gammie 2001). This may 
provide a local source of energy and maintain a Toomre pa- 
rameter Q t ~ 1 (Toomre 1964). In a marginally gravitationally 
unstable thin disk, the non-axisymmetric perturbations lead to 
angular momentum transport that can be described by an ef- 
fective viscosity (i.e., Paczynski 1978; Shlosman et al. 1989, 
1990; Rafikov 2009). Self-consistent models of geometrically 
thick disks in AGNs at parsec scales call for the effects of self- 
gravity and radiation physics in three dimensions. This difficult 
endeavor can be the subject of future research. Another compli- 
cation is that the problem of parsec- scale winds may be intrin- 
sically connected with the problem of accretion at parsec scale. 
Relaxing the assumption of an equatorial thin disk and treating 
accretion adopting an effective a prescription is a possible next 
step in the development of the current model. 

The assumption of a razor-thin disk as a source of matter 
underpins the current study. Certainly such a distribution of 
matter is in many respects most unfavorable for the formation 
of radiation-driven flow, owing to the geometry of X-ray 
illumination. Thus, we expect that if the wind does form for 
certain T in the thin-disk case, then we may expect that it 
will also form in less idealized circumstances. To relax the 
assumption of a razor-thin disk, one would need to abandon the 
assumption of equatorial symmetry and to include an accretion 
disk in self-consistent consideration. 

Starting from Paper I, we see that much of the massive wind 
does not leave the potential well of the BH. This conclusion is 
confirmed through Paper II, as well as in the present work. The 
precise fate of this gas is beyond the scope of these works, but 
the results suggest the following: if the accreted gas is rejected 
by the BH, at least some of it will return for additional attempts. 

We conclude with a mix of results and ideas driven by our 
current studies: 

1 . The distribution of radiation force inside the dusty compo- 
nent of the flow depends on the shape of the photosphere 
where X-rays are converted into IR and can only be calcu- 
lated in multi-dimensional simulations. This work includes 
both X-rays and IR in 2.5D time-dependent, hydrodynamic 
simulations. 

2. We do not observe a rotationally supported quasi-static 
torus in our simulations, although our simulated winds are 


likely to resemble some of the observed properties of such 
tori. 

3. The flow in the current simulations is more complex and 
time-dependent than that of our previous studies, where we 
did not include X-rays explicitly. 

4. If geometrically thick dusty obscuration develops, then a 
hot photoionized flow with velocities of 100-700 km s _1 
accompanies it. X-ray warm absorbers are the evaporative 
flow both originating from the disk and also evaporated 
from such a cold, dusty component. 

5. We speculate that large-scale motions with v z = ±(afew)x 
100 km s -1 originating from a spatially varying wind such 
as we have calculated may be seen in maser emission in 
nearby bright type II AGNs. 
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APPENDIX A 

Frequency-integrated moments E , F, which appear in the 
above set of Equations (l)-(4), are obtained by calculating an- 
gular moments from the frequency-integrated specific intensity, 
/( r, £2, v, t ): 

E(r,t) = -J dv <j) d£l I(r, £1, v, t), (Al) 


F(r, t) 



d£2nl(r, £2, v, t). 


(A2) 


The frequency-independent radiation pressure tensor P is found 
from 


P(r, t) 


1 

c 



d£liml(r , £2, v, t). 


(A3) 


APPENDIX B 


Dust directly reprocesses X-rays to IR, and we approximately 
take this into account. The amount of energy absorbed by dust 
in a volume d 3 x is: 

1 dE v C 

d 3 x~dt L=and / °Jv*d£2dv ~ an d a x cu x , (Bl) 

where u x is the energy density of X-rays, a x = Xx/^d, and the 
number density of dust grains reads n d = fdP/wid, where m d 
is the mass of a grain and/d is the dust-to-gas mass ratio. We 
make the further simplifying assumption that the dust grain is 
being instantaneously heated to the effective temperature T dx . 

The contribution from the dust to the IR energy density 
(or, equivalently, to the temperature T r ) is calculated from the 
following equation: 


= n d c a d (r d>x - T t ) . (B2) 

Assuming that 7d )X is constant over the time interval [t, t + dt], 
the update for radiation is Td, x (t + dt) = 7d ;X (l — e ~ dt/td ) + 
— e ~ dt / td Td, x (t). Note that E = aT r 4 . Equation (B) is used in an 
operator- split fashion to update E\ prior to this all other terms 
in Equation (4) are updated. 
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